Multiple channelling single-photon emission with scattering holography designed metasurfaces

Channelling single-photon emission in multiple well-defined directions and simultaneously controlling its polarization characteristics is highly desirable for numerous quantum technology applications. We show that this can be achieved by using quantum emitters (QEs) nonradiatively coupled to surface plasmon polaritons (SPPs), which are scattered into outgoing free-propagating waves by appropriately designed metasurfaces. The QE-coupled metasurface design is based on the scattering holography approach with radially diverging SPPs as reference waves. Using holographic metasurfaces fabricated around nanodiamonds with single Ge vacancy centres, we experimentally demonstrate on-chip integrated efficient generation of two well-collimated single-photon beams propagating along different 15° off-normal directions with orthogonal linear polarizations.

Holography is an approach used for recording and reconstructing arbitrary wavefronts.A hologram is produced by recording an interference pattern between a reference wave and a signal wave, typically by using a photosensitive thin film.When illuminated with the reference wave, the hologram, representing usually an amplitude or phase-encoded transparency, reproduces (with certain accuracy) the signal wave phase and amplitude distribution.The intensity interference pattern can also be calculated and thereby imprinted into a thin film (varying either its amplitude transmission or thickness), resulting in a computer-generated hologram.In this section, we describe our holographic method, scattering holography S1 , developed for designing metasurfaces that produce far-field beams of predefined directions and polarizations.
We discuss in the following the case of a cylindrically diverging SPP wave, which is assumed to be excited by a quantum emitter (QE), used as a reference wave interacting with a (computergenerated) metasurface and scattering into two orthogonally linearly polarized beams separated in their direction of propagation.

S1.1 Holographic approach: shaping QE emission
First, we present a general method of recording the holographic pattern and reconstruction of the signal wave using SPP as a reference wave.Then, we elaborate on the case of QE-excited cylindrically diverging SPP and provide corrections to the model.As a starting point, we consider the interference between the reference SPP wave E spp and the signal wave E s at xy-plane [z=0, Figure S1(a)].The interference pattern can be written as: where E 0 spp and E 0 s denote the complex vectorial electric field amplitudes, while k s and k spp represent wavevectors for the signal and reference SPP waves, respectively.All quantities in equation S1 are functions of coordinates (x, y).The in-plane radius vector is denoted by r xy .Here, the in-plane electric foeld component of QE-excited SPP E 0 spp is oriented radially: E 0 spp (x, y) = E 0 spp (x, y) (cos(ϕ), si n(ϕ)), where ϕ is the angle between x-axis and the SPP propagation direction k spp .The in-plane components of vectors k spp and r xy are parallel to each other, therefore k spp • r xy = k spp r x y .
We consider the hologram to be composed of isotropic, non-interacting dipolar nanoscatterers, for example, spherical nanoparticles S1 .The volumes of these nanoparticles, and hence their scalar polarizabilities, are considered to be proportional to the intensity of the interference pattern: V sc (x, y) ∼ I (x, y).The size of scatterers is assumed to be much smaller than the operating light wavelength so that one could employ the electric-dipole approximation when considering radiation scattering by the hologram nanoparticles.
When such a hologram of distributed nanoparticles is illuminated with the reference SPP wave, the scattered field in the immediate vicinity of a nanoparticle is proportional to the incident field and the nanoparticle polarizability, which is proportional to the nanoparticle volume and thus to the local intensity of the interference pattern: The first term of the scattered field represents the evanescent field with k spp > k 0 which doesn't propagate away from the hologram's xy plane.The second term of E sc is also a characteristic evanescent wave because the associated wavevector magnitude is larger than wavenumber k 0 of a free propagating wave, 2k spp − k s > k spp > k 0 .As a result, the only reconstructed wave, which propagates away from the hologram surface, is the reconstructed signal wave represented by the last term: E sc .This propagating away scattered field propagates in the direction of the signal wave [Figure S1(b)], but has the polarization of the reference wave:

S1.2 Holographic approach: channeling QE emission into two beams
In this section, we consider the case of holographic reconstruction of two beams with orthogonal polarizations.
V The interference between the reference E spp and the two signal waves E s1 , E s2 [Figure S2(a)] can be calculated as follows: The field reconstructed by the holographic array of nano-scatterers can then be written down: The radiative part of the reconstructed field distribution, E sc (x, y), contains two terms (terms with k s1 , k s2 ) representing the signal waves [Figure S2(b)].Additionally, a part of the radiated light is diffusely scattered, also propagating away from the surface (the last two terms).However, in our case of orthogonally polarized signal waves, the dot products of the two signal field amplitudes become zero.Therefore, once again only the signal waves are propagating away in the process of hologram reconstruction: S-7

S1.3 Cylindrically diverging SPP as a reference wave
Now, we elaborate further on the form of the reference SPP wave, which is a QE-excited radially divergent SPP wave: where α = 1 2L spp , and L spp is the SPP propagation length.The SPP wavenumber where N spp is SPP effective mode index.
One can observe that the SPP attenuates during its propagation due to its radial divergence ( 1 r x y ) and absorption (exp(−αr x y )).In the process of the hologram reconstruction, the SPP is also depleted owing to its out-of-plane scattering by hologram nanoparticles.The above factors limit the domain of interference and even more so the hologram area efficiently used in the reconstruction.Considering the above SPP field (eq.S7) used as the reference wave during the hologram recording and reconstruction, the reconstructed signal wave term (similar to eq.S3) becomes: Here, the inverse radially and exponentially decaying factors appear, limiting the hologram spatial domain used for the signal reconstruction to a small area around the QE and making faithful reconstruction rather problematic.This problem can however be circumvented, if we compensate these decaying terms at the stage of hologram recording, when the latter is implemented by calculating the interference pattern and fabricating the corresponding (computer-generated) hologram, by using the artificially growing SPP field, i.e., calculating the interfence pattern , where: This compensation contrivance allows one to reconstruct the signal wave, which is relatively S-8 faithful as far as the amplitude and phase distribution near the array of hologram nanoparticles is concerned (although leaving the problem of polarization reconstruction unresolved): where the SPP decay factors are compensated, and the spatial domain of the reconstruction is not Our holographic approach developed for designing the QE-coupled metasurfaces allows one to generate the QE emission with arbitrary wavefronts.We have used this approach to calculate the metasurface pattern for generating two angular-resolved emission channels with two orthogonal linear polarizations.First, the xy-plane field components are calculated for the desired signal beams.Second, the interference pattern is calculated as described above.Finally, the obtained interference distribution I (x, y) is discretized into an array of scatterers, which constitute the dielectric metasurface.

S2 Calculation of metasurface pattern for reconstruction of two signal waves of orthogonal linear polarizations
As described in the previous section, to calculate the metasurface pattern, we first calculate the interference pattern of the two signal waves and the artificial SPP wave and further discretize it into a binary pattern with pixels either filled with dielectric or kept empty.The wavelength of the light is set to λ 0 = 602 nm to match the GeV emission peak.

S-9
The calculation is two-dimensional, and all the fields are projected onto the xy-plane.The field components, which are normal to the hologram (metasurface) plane, are disregarded, because their scattering along directions close to the normal is negligibly small (within our electric-dipole approximation).The obtained metasurface pattern is also two dimensional being associated with the two-dimensional intensity interference of signal and reference waves.
The considered system is a planar interface between silver coated with 30 nm SiO 2 and air oriented in the xy-plane.
Two Gaussian beams of orthogonal linear polarizations are considered as signal waves: The first signal beam has the polarization along the x-axis and the propagation direction k x s in the zy-plane, beam waist ω 0 .The deflection angle from the z-axis is θ = 15 • as illustrated in As discussed in the section S1.3, artificially increasing the SPP wave(eq.S9) is used to compensate for the SPP divergence and absorption (and eventually scattering) and, in doing so, to eliminate the spatial limitation of the domain of its interference with the signal waves.The artificially growing SPP wave is radially divergent from the origin, (x, y) = (0, 0): where, k spp = k 0 N e f f is SPP propagation constant, k 0 = 2π/λ 0 the wave vector of free-space propagating wave, λ 0 = 602 nm.r xy denotes the radius vector.The SPP effective mode index N e f f is calculated as follows: The effective mode indices for the Ag/SiO 2 /HSQ and Ag/SiO 2 /air interfaces are numerically calculated separately as N e f f (ai r ) = 1.12 and N e f f (H SQ) = 1.52, respectively, resulting in N e f f = 1.32.The filling factor χ was fixed at 50%, denoting the dielectric/air ratio of the optical constants of silver were performed for a comparison (Figure S6).
Next, the interference between the reference E a,spp and the two signal waves E x s , E y s is calculated as described above.Note that intensity interference patterns are typically composed S-11 of interference fringes.According to our scattering holography approach for recording the hologram, one should therefore arrange (spherical) nanoparticles along the lines (interference fringes) varying their volumes in accordance with the local intensity levels.Such an approach although being feasible (and in line with the theoretical foundation described above) is extremely demanding as far as practical implementation is concerned.Targeting the conventional fabrication, for example with the electron-beam lithography, we consider the nanoparticles being fused into continuous nanoridges of the constant height but with the local width being proportional to the local intensity in the interference pattern.This approach is adopted throughout the present work, resulting in high quality and efficient reconstruction of signal waves by QE-excited SPP being scattered out with the designed metasurfaces, as demonstrated in this work. S-12

S2.1 Initial design
For the first approximation, the hologram recording and reconstruction processes can be built under the following assumptions: 1. Constant SPP amplitude, neglecting the SPP field attenuation due to divergence factor 1 r x y , absorption [exp(−αr x y )] and depletion owing to its scattering by hologram nanoparticles.
The interference patterns I (x, y) that emerge from the interference between E a0,spp and projections of signal wave fields E x s,pl ane and E y s,pl ane are shown in Figures S4(a-b).The final interference pattern of all three waves is shown in Figure S4(c).
Next, we discretize the obtained hologram pattern to generate the design pattern suitable for the fabrication with the electron-beam lithography (as discussed in the preceding section).The field intensity distribution is normalized and lies within the rang e [0-1].The threshold value is defined by the structure filling factor (in our case χ = 0.5).In this case, half of the metasurface area is filled with material.Therefore, half of the area of the interference pattern should be below the threshold and half should be above.Thus we find the value of threshold and discretize the structure treating all pixels larger than threshold as "1" and, oppositely, all pixels smaller than threshold as "0."The binary image is shown in Figure S4(h).All "1" pixels correspond to the positions of noninteracting isotropic scatterers.In our case, these points were filled with the HSQ.The calculated far-field distributions presented in Figure S4(e-g) show the generation of two separate beams with orthogonal polarizations.

S2.2 Optimized design
To elevate the beam quality in terms of collimation and polarization clarity, the Gaussian profile of the signal beams and depletion of the reference SPP wave are considered.
The Gaussian profile of the signal waves is represented by the exp(−r 2 /ω 2 0 ) term in field E x(y) of the beam (eq.S14).The beam waist, ω 0 , becomes an important optimization parameter because it limits the spatial domain of the interference pattern calculation, as shown by Kan et al.S1 .
As discussed above, to compensate and enlarge the interference calculation domain, we utilize an artificial SPP field that grows with distance from the QE [S15]: This procedure results in improved beams collimation, as shown in Figures S6(d-f).Comparison of the beam profiles clearly shows higher peak intensity and better collimation of the optimized beam (c).Table (d) shows results of optimization in terms of primary metasurface parameters.The comparison of the Collection Efficiency (CE) for 0.3, 0.5 and 0.9 NA objectives is displayed in table (e).

S2.3 Performance comparison of initial and optimized metasurface design
A comparison of the emission characteristics of the initial and optimized metasurface design is presented in Figure S6.The optimized structure is designed with consideration to Gaussian profiles of the signal beams and compensation for SPP decay.As a result, the reconstructed signal beams are significantly better collimated, featuring lower overlap and stronger peak intensities The first parameter of comparison is the quantum efficiency calculated (for details of calcula-tion -see section S3) using the Johnson & Christy S2 and Palik S3 optical constants for the silver layer of the substrate.The optimized structure exhibits 63% and 45% quantum efficiencies for Johnson & Christy and Palik data, respectively.The quantum efficiency value is lower for the Palik case because of the higher absorption in silver.As expected, the external quantum efficiency of the optimized device is not enhanced because there is no change in the energy contained in the SPP wave, and the energy is only redirected in a more collimated way.The next parameter is the overlap integral, which is a figure-of-merit for beams' separation.The lower the value of the overlap integral, the better the splitting of the beams.The overlap integral is calculated as where I x and I y are the intensities of the x-and y-polarized emission, respectively, and integration is performed over a 30 • emission cone, which is projected to the far-field as a disc.
The integration area contains both emission lobes.The peak area measurement allows us to conveniently compare the beam collimation efficiency.The area of the peaks is calculated by counting all the pixels of the far-field emission distribution with an intensity higher than the 10% threshold.The divergence angle θ X (Y ) 1/2 is calculated as angular y-polarized beam full width at half maximum along the x(y) direction.We consider only the y-polarized beam because beams are symmetrical relative to the structure diagonal.
From the comparison, we observe that the beams divergence is improved from 4.3 . The overlap integral is reduced from 17.4% to 12.8%.This leads to a significant enhancement in collection efficiency for low-NA objectives (0.3, 0.5 NA) (S6(e)).The CE is calculated as the ratio of power emitted in a certain NA to the power emitted in 90 • (NA = 1).

S3 Quantum efficiency calculation
For a quantum emitter coupled to a metasurface, the quantum efficiency can be written as: where η QE denotes quantum efficiency of the coupled system, η IQE denotes the internal quantum efficiency of the quantum emitter, and η EQE denotes the external quantum efficiency of the coupled system.η IQE depends on the properties of the emitter and can vary between 0 and 1 for different emitters.EQE depends on the device, in our case, a metasurface.With η IQE assumed unity for an emitter, the quantum efficiency of the coupled system will be equal to external quantum efficiency.In this section, we present a method that we have utilized to calculate the external quantum efficiency of our device.The external quantum efficiency was calculated numerically using an FDTD software.η QE was calulated as the ratio between the power radiated into outgoing waves, propagating within a 64 • cone (collected by an NA = 0.9 objective) and total power emitted by QE in the presence of the plasmonic environment (metasurface).
The power generated by the dipole in the presence of the plasmonic environment was S-20 calculated by integrating the Poynting vector over a 100x100x30 nm (x,y,z) box enclosing the point dipole.The power outcoupled by the metasurface to the propagating waves was calculated by integrating the Poynting vector over the cuboid cap, covering the metasurface, and placed at a reasonable distance (300 nm) from the top of the metasurface to avoid the influence of the near-field.The size of the cap was 20×20×0.35µm, covering the metasurface with a diameter of 17 µm.To be consistent with the experimental measurements, the size of the cap was adjusted to count only the radiation emitted in the NA=0.9 (64 deg).This procedure resulted in η QE = 0.63 external quantum efficiency.Predictions obtained using a similar method were experimentally verified before S4,S5 .It should be noted that the efficiency can be further increased by using high-refractive index metasurface ridges to boost their out-of-plane scattering and high-quality monocrystalline metal films to reduce the ohmic losses S6 .Thus, simply by changing the ridge material to titanium dioxide (TiO2), the external quantum efficiency is expected to exceed 80%, even without further optimization of ridge parameters as described for the bullseye pattern S5 .
The quantum efficiency of the coupled system depends on the scattering and propagation loss of SPPs.Below, we present a method to separately estimate the effects of scattering and propagation loss on quantum efficiency.
The SPP energy is either scattered by metasurface ridges or absorbed by metal because of the ohmic losses: where P 0 is part of the dipole emission power that is coupled to the SPP mode.
So, the outcoupling efficiency can be defined as: To estimate P sc at , we need to consider the process of the SPP scattering by the metasurface ridges depending on the distance from the SPP source.We consider 2D metal-dielectric interface S-21 (xy-plane).The field of the radially divergent SPP propagating on the clean metal-dielectric interface (without metasurface) is described by: where α abs = 1 2L spp , with L spp being the SPP propagation length (over which the SPP intensity decreases by 1/e), and k spp is SPP wavevector with corresponding wavenumber k spp = 2π λ N spp , with N spp being the SPP effective index at the operating wavelength λ.The radial divergence is represented by the term 1 r , where r is the radius.Here α abs is defined for the ohmic losses in metal.
Therefore, the SPP intensity is: In the presence of the metasurface the SPP field is also depleted due to scattering by the dielectric ridges.Therefore, the attenuation coefficient of the SPP interacting with metasurface ridges should be changed to: and the SPP intensity will decay faster: Assuming that power scattered in the thin ring with the thickness d r is proportional to the SPP intensity: The total scattered power can be calculated by integrating d P sc at over the metasurface area: The result is: The dependance P sc at (L) can be obtained by estimation of the outcoupled power for metasurfaces with increasing diameters in numerical simulations.Using this function to fit the simulated data allows us to estimate the decay coefficient α t ot = 0.13 .
Following a similar procedure, the part of the energy absorbed in metal can be written as: The SPP decay coefficient α abs = 0.03 was calculated analytically: In our case λ = 602 nm and N spp = 1.28+i •2.77•10 −3 calculated for Ag-SiO 2 -HSQ-air interface S-23 with HSQ filling factor 0.5, which corresponds to the metasurface presented in the manuscript.
The SPP decay coefficient attributed to scattering is calculated as: This results in 78% outcoupling efficiency (α sc at = 0.1, α t ot = 0.13): ε sc at = P sc at P 0 = P sc at P sc at + P abs = α sc at α sc at +α abs = α sc at α t ot = 0.78 (S33) The external quantum efficiency is a product of efficiency of the dipole coupling to SPP η spp , scattering of SPPS by the metasurface ridges η sc at and collection efficiency of the objective η C E , that is, First, the excitation efficiency of SPP waves is η spp > 0.8 for proper orientation of the dipole placed at a reasonable distance from the metal surface S7 .This value can also be derived using an analytical model as well S8 .Second, the scatterd part of the propagating SPP, by the dielectric ridges of the metasurface is calculated to be η sc at ≈ 0.78, the remaining 22% are absorbed by the metal.The collection efficiency, defined as the fraction of the total outcoupled power collected by the objective of NA=0.9, was numerically calculated to be η C E = 0.98.Taking all the aforementioned factors into account, the calculated total external quantum efficiency of is: η QE = 0.8 * 0.78 * 0.98 = 0.61, which is close to the value simulated numerically.
Using high-refractive index metasurface ridges, for example made of titanium dioxide, would further increase the contribution of the scattering SPP loss, making the SPP absorption loss contribution negligible and reaching the total external quantum efficiency of 0.8 (as was obtained for the bullseye pattern S5 .We fabricated our samples in following five steps:

S5 Metasurface fabrication
1. Substrate preparation: A polished Si wafer is coated by thermal evaporation of 3 nm Ti followed by 150 nm Ag, and another 3 nm Ti.Radio frequency magnetron sputtering of 30 nm SiO 2 layer is done to protect silver from oxidation.Ti layer is used for better S-26 adhesion between Si and SiO 2 to metal.Deposition rates for Ti/Ag/Ti/SiO 2 structure are 0.1/1/0.1/1.5 Å/s, respectively.Chamber pressure used for depositions is ∼ 5 • 10 −6 mbar.The positioning accuracy is ≈50 nm.The positioning procedure is described in previous work S4 .The procedure of GeV-ND fabrication is described elsewhere S9 ).

Figure S1 :
Figure S1: Concept of scattering holography with cylindrically diverging SPP generated by a quantum emitter (QE).Calculation of the interference pattern intensity generated by a signal wave E s and SPP reference wave E spp in the hologram xy-plane (a), and reconstruction of the signal wave by scattering of the reference SPP wave by nanoparticle distribution matching the interference pattern (b).

Figure S2 :
Figure S2: A schematic of vectorial scattering holography with cylindrically diverging SPP generated by quantum emitter (QE).Calculation of the interference pattern intensity generated by two signal waves E s1 , E s2 and SPP reference wave E spp in the hologram xy-plane (a), and reconstruction of the signal wave by scattering of the reference SPP wave by the nanoparticle distribution matching the interference pattern (b).
explicitly limited.The SPP depletion because of the scattering by the hologram nanoparticles can be compensated by inserting the additional decay term α sc at in the attenuation coefficient α * = α+α scat .It should however be borne in mind that the reconstruction domain is implicitly limited by the requirement of hologram nanoparticles being of subwavelength sizes: the exponential growth of the SPP field used in the hologram recording requires a similar growth of the hologram nanoparticle volumes.

Figure S3 .
Figure S3.The second signal beam has the polarization along the y-axis and the propagation direction k y s in the zx-plane with a 15 • deflection from the z-axis.The projections of the signal waves in the xy-plane are as follows:

Figure S3 :
Figure S3: Projection of the X-polarized signal wave onto the XY plane.For interference pattern calculation the signal wave ( field E 0 polarized along X-axis (blue), wavevector k x 0 ) is projected onto XY plane.The Y-polarized signal wave is projected in the similar manner.Angle θ defines the inclination of the signal wave in respect to the normal to the hologram plane and direction of the reconstructed wave.
metasurface.The fixed filling factor allowed simplification of the calculations by utilizing the a constant SPP mode index for the structure.The refractive indices of the dielectrics were set to n(SiO 2 ) = 1.45 and n(H SQ) = 1.41.The SPP decay factor α abs = 1 2L spp = I m{k spp }.In our experiment, the SPP propagation length is L spp =17 µm and α abs = 0.03.The Johnson&Christy S2 silver optical constants were used for pattern generation.Also, calculations based on Palik S3

Figure S4 :
Figure S4: Initial design of the hologram metasurface.Interference between SPP and X-or Y-polarized beam fields on the sample surface is shown in (a) and (b), respectively.The final interference pattern (c) between SPP, X and Y polarized beams.(d) The cross section of (c) with the threshold used for discretization.(h) Binary image of metasurface suitable for the electron beam lithography.(e-g) Simulated far-field patterns representing emission distributions inside the 30 • cone.

Figure S5 :
Figure S5: Optimized design of the hologram metasurface.(a) The interference pattern between artificial SPP, x-and y-polarized Gaussian beams.(b) The cross section of (a) with the threshold used for discretization.(c) Binary image of metasurface suitable for the electron beam lithography.(d-f) Simulated far-field patterns representing emission distributions inside the 30 • cone.

Figure S6 :
Figure S6: Comparison of the performance of the initial and optimized metasurfaces (a) Design of the metasurfaces.(b) Improvement of the beam collimation in resulting far-field emission patterns.Comparison of the beam profiles clearly shows higher peak intensity and better collimation of the optimized beam (c).Table (d)  shows results of optimization in terms of primary metasurface parameters.The comparison of the Collection Efficiency (CE) for 0.3, 0.5 and 0.9 NA objectives is displayed in table (e).

Figure S7 :
Figure S7: Metasurface size comparison.(a) Top view of metasurfaces with corresponding diameters of 5, 10, 15, and 17 µm.(b) Calculated far-field patterns representing emission distribution inside 30 • cone.The increase of the structure diameter enhances the beam collimation.(c) Comparison table of metasurface performance.

Figure S8 :
Figure S8: Calculation of external quantum efficiency.The power radiated by the point dipole was calculated by integrating the Poynting vector over the 100x100x30 nm (x,y,z) box around the dipole (a).The power of the outcoupled emission was calculated by integrating the Poynting vector over a 20x20x0.35µm cap on top of the metasurface (b).The ratio of the two powers is the external quantum efficiency.The thickness of the SiO 2 layer is 30 nm.The thickness of the Ag layer is 320 nm.

Figure S9 :
Figure S9: Calculation of the scattered SPP power by dividing the metasurface domain into rings with thickness d r .Each ring located at a distance r from the SPP source (point dipole) scatters power d P sc at ∼ α sc at I spp 2πr d r .Integration of the d P scat over the metasurface domain r = [0 L] results in total scattered power by metasurface P scat .

2.
Gold markers fabrication by electron beam lithography: PMMA A2 positive resist is spincoated at 1580 rpm for 45 sec and baked at 180 o C on a hot plate for 2 min to form a ≈100 nm layer.10x10 µm 2 crosses are patterned using an electron beam lithography (EBL) system and developed for 35 sec in 1:3 MIBK:IPA (Methyl isobutyl ketone:Isopropyl alcohol) and rinsed in IPA for 60 sec.Ti(3nm)/Au(35nm) layers are thermally deposited in a vacuum chamber with deposition rates 0.1/1 Å/s and pressure 5•10 −6 mbar.The unexposed PMMA is dissolved in acetone during a 12 hours lift-off procedure.Samples are ultrasonicated in acetone for 1 min and washed with IPA and DI water.Gold is used as a material for markers because it is chemically stable and provides sufficient contrast for optical microscopy and fluorescent mapping.3.Nanodiamond deposition:To stabilize nanodiamonds on the surface we cover the substrate with an adhesive Poly(allylamine hydrochloride).Further, a water dispersion of Nanodiamonds containing Germanium vacancy centers (GeV-NDs) with average size ∼100 nm size is spin-coated on the sample.The positions of diamonds are determined relative to the gold markers using dark-field microscope images shown in FigureS11(a).

4.
Nanodiamond choice: The crucial part of QE-SPP coupled metasurface fabrication is to ensure that the QE radiative transition dipole is predominantly oriented perpendicular to the sample surface, since this results in a high QE-SPP coupling efficiency along with axially symmetric propagating SPPs S4,S5,S10 .To identify NDs with properly oriented GeV centers, we used a strongly focused and radially polarized excitation beam, producing a strong E z component of the electric field in the focal point and consequently efficiently exciting GeV centers having a predominant out of plane component of the dipole moment.In the S-27 last step, we fabricate HSQ holographic metasurfaces around individual GeV-NDs using their coordinates determined relative to the markers to properly conduct electron-beam lithography (EBL). 5. Hydrogen silsesquioxane (HSQ) structure fabrication: Negative resist HSQ (Applied Quantum Materials, Canada, 8% solution in MIBK) is spin-coated at 1200 RPM 60 sec and baked at 160 o C on a hot plate for 2 min to form a ≈200 nm layer.The metasurface structures are patterned around GeV-NDs by EBL and developed in 25% TMAH (Tetramethylammonium hydroxide) for 4 min.Samples are rinsed for 1 min with DI water.Images of a fabricated structure are shown in Figure S11(b) and Figure S12.